A Critical Look at Cosmological Perturbation Theory Techniques 
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Recently a number of analytic prescriptions for computing the non-linear matter power spectrum 
have appeared in the literature. These typically involve resummation or closure prescriptions which 
do not have a rigorous error control, thus they must be compared with numerical simulations to 
assess their range of validity. We present a direct side-by-side comparison of several of these analytic 
approaches, using a suite of high-resolution N-body simulations as a reference, and discuss some 
general trends. All of the analytic results correctly predict the behavior of the power spectrum at 
the onset of non-linearity, and improve upon a pure linear theory description at very large scales. 
All of these theories fail at sufficiently small scales. At low redshift the dynamic range in scale where 
perturbation theory is both relevant and reliable can be quite small. We also compute for the first 
time the 2-loop contribution to standard perturbation theory for CDM models, finding improved 
agreement with simulations at large redshift. At low redshifts however the 2-loop term is larger 
than the 1-loop term on quasi-linear scales, indicating a breakdown of the perturbation expansion. 
Finally, we comment on possible implications of our results for future studies. 
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I. INTRODUCTION 



The character and evolution of the large-scale struc- 
ture of the Universe has been the subject of much re- 
search in recent decades. As it is currently understood, 
large-scale structure grows through a process of gravi- 
tational instability starting from a nearly scale-invariant 
spectrum of Gaussian fluctuations at early times. On 
very large scales the matter distribution of our universe 
today is well modeled by linear perturbation theory. On 
scales below about lOMpc, on the other hand, the dy- 
namics are highly non-linear and we must resort to direct 
numerical simulations of the N-body problem to under- 
stand the clustering of matter or its tracers. 

On intermediate, or quasi-linear, scales there is the 
possibility that the matter distribution may be mod- 
eled analytically by extending perturbation theory be- 
yond linear order. This possibility has received renewed 
attention recently due to the interest in using baryon 
acoustic oscillations as a probe of the expansion history 
of the Universe and of the nature of dark energy [l| . Since 
the baryonic features are at large scales (0(100) Mpc) it 
is plausible that higher order perturbation theory could 
model subtle corrections to the linear result with some 
accuracy. More generally, investigation of perturbation 
theory may allow some improvement in theoretical pre- 
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dictions for the next generation of very large surveys. 

Consequently, a number of new ideas have been in- 
troduced in recent years for computing statistical prop- 
erties of the matter distribution, most importantly the 
2-point function or power spectrum. Regrettably these 
approaches involve uncontrolled approximations, provid- 
ing no simple way of estimating the theoretical uncer- 
tainty. Since perturbation theory is expected to fail on 
sufficiently small scales, the domain of validity of any 
particular approach is therefore unclear, and the only 
known way to test their accuracy is to compare their 
predictions with the results of N-body simulations. In 
the past this has been done on a case-by-case basis, with 
one theory tested for one cosmology against one suite 
of N-body simulations, focusing primarily on the power 
spectrum. Recently there have been some attempts to 
compare multiple theories simultaneously 0, Q j or to ex- 
amine statistics other than the power spectrum 0, S @1 • 
However a comprehensive comparison has been lacking, 
and with the recent proliferation of analytic techniques it 
is natural to ask how well these theories actually perform. 
With near-future observations potentially depending on 
these techniques and with recent advances in N-body al- 
gorithms and computing power, it is timely to revisit this 
issue. 

In this paper we present a direct comparison of sev- 
eral recent analytic predictions for the clustering of mat- 
ter on quasi-linear scales. We restrict our attention to 
the matter fluctuations, because very few of the existing 
treatments can handle biased tracers such as dark matter 
halos and galaxies. We use modern, high-resolution N- 
body simulations as our reference points, which provide 



2 



highly accurate @, M, (though computationally expen- 
sive) estimates for statistical observables of the matter 
distribution. By comparing the analytic predictions for 
two cosmologies, one close to the current best-fit model 
and one more extreme, we are able to judge the relative 
merits of each approach. 

The paper is organized as follows. In Section |TT] we 
start by reviewing the dynamical equations that govern 
the evolution of the matter distribution and discuss the 
relevant statistical quantities that one may compute. We 
then continue by summarizing the different analytic ap- 
proaches we consider in this paper. In Section [Ml we de- 
scribe the N-body simulations that are used as a reference 
point for the comparison. In Section Hvl we plot the var- 
ious approaches together, discuss qualitatively how well 
they agree with simulations, and propose several ways to 
quantify this agreement. We discuss the results of this 
comparison in SectionfVl and make some closing remarks 
in Section IVT1 



II. ANALYTIC METHODS 

We start by reviewing the different analytic meth- 
ods we consider - our goal is not to provide a compre- 
hensive description of each method, but to provide an 
overview and highlight the relationships between the dif- 
ferent methods. 



A. Dynamics and Linear Theory 

By far the most popular approach to an analytic de- 
scription of large-scale structure is to approximate the 
matter distribution as an irrotational fluid, characterized 
by a density constrast 8(x) = p(x)/p — 1 and a peculiar 
velocity divergence 6{x) = V -v(x). The fluid equations, 
in Fourier space, are then (see Appendix lAl for a detailed 
derivation) , 
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Here dr = dt/a(t) is conformal time, f2 TO (r) = 
P~( T ) / Pcrit(i~) , and TL = aH is the conformal Hubble pa- 
rameter. (Note that we adopt the Fourier transform con- 
vention that puts the (2ir) 3 in the wavevector integral. 
We also omit the tilde that is usually used to decorate 
Fourier space quantities.) The non- linear nature of these 
equations is manifest in the mode-coupling integrals. 
Working to linear order in 6 and 9, we obtain 

S L (k;z) = ^-5 i (k) (2) 

and 

9 L (k;z) = -H(z)f(z)^-8 i (k), (3) 

where Si is the density contrast at some early time 
Zi when linear theory is certainly valid, D is the lin- 
ear growth function (normalized to 1 at z = 0), and 
/ = din D/dhia. At early times Q m w 1 and D oc a. 
For convenience we define Jo to be the linear density 
contrast today, i.e. 8o(k) = &L{k\z = 0). When con- 
venient we also follow common practice and use rj = In D 
as a time variable; for brevity we often suppress the time 
dependence of quantities altogether. It is further con- 
venient to group S and 9 into a 2-componcnt vector, 
$„(fe) = (5(k), -9(k)/Hf) which is proportional to (1, 1) 
in linear theory. 



B. Statistical observables 

Inflation predicts, and observations have confirmed, 
that the initial fluctuations are predominantly adiabatic 
[lol ] , almost scale- invariant , and very close to Gaus- 
sian [ll|. Under the assumption that the initial field 
is Gaussian all expectation values of moments of the 
evolved density and velocity fields can be expressed as 
integrals over the linear theory power spectrum. For ex- 
ample, the evolved 2-point function 

(2ir) 3 5 D (k + k')P ab {k) = ($ a (fe)$ b (fc')> , (4) 

whose components are all equal to Pl (k) in linear theory, 
can be expressed as integrals over n powers of Pl in n th 
order perturbation theory (e.g. Eq. (^201) ). 

In general, to give a complete statistical description of 
the matter distribution at a given time, one would need 
to specify the entire hierarchy of connected n-point cor- 
relators. For initially Gaussian fields which are close to 
linear, the higher order connection functions are small 
and have been compared to simulations in [l2j • We shall 
confine our attention to the 2-point function in this pa- 
per. 

The non-linear propag ator ([H E3; 

also known as 

the response function [J|) measures the correlation be- 
tween the evolved field <f> a (fc; rj) and the initial conditions 
^a(k;rji). It is formally defined as a functional deriva- 
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tive, 

though its significance is easier to understand from the 
relation 

($ Q (fc; ri)$b(k'; m)) = G ac (k; 77, rji) ($ c (fc; rn)$ b (k'; %)) , 

(6) 

which we shall take as a definition henceforth. At early 
times or at large scales there is near-perfect correlation 
{Gab ~ 1), but Gab —* on small scales as non-linear 
evolution washes out the initial conditions. 

Because we will make reference to it later, we also in- 
troduce here the quantity 

1 f°° 

which characterizes the scale at which non-linearities be- 
come important. In the Lagrangian formalism (see be- 
low) S 2 gives the variance of each component of the linear 
(or Zel'dovich) displacement. 

C. Beyond Linear Theory 

The program is now to compute the statistics of the 
evolved density field in terms of the initial density field. 
This is simple in principle but difficult in practice, be- 
cause the equations of motion are both non-linear and 
non-local (in both configuration space and Fourier space) . 
Non-linearity forces one to seek a perturbative solution, 
since exact solutions to Eqs. ([T]) (even if they could be 
found) could not be combined to construct a realistic 
solution. A straightforward perturbative approach is 
hampered by computational costs, as non-locality implies 
that higher order terms involve mode-coupling integrals 
of ever higher dimension. 

This situation has prompted a study of higher-order 
methods for statistical observables like the power spec- 
trum. Many of these methods were borrowed from other 
areas of physics (notably particle physics and fluid me- 
chanics [l3) where they achieved mixed success. We re- 
view these below, highlighting the relationships between 
the different methods; the methods we consider are sum- 
marized in Table HT1 

The most straightforward approach is to define a series 
solution to the fluid equations in powers of the initial den- 
sity field Si (or cquivalently, the linearly evolved density 
field, So). This is the basis behind standard perturba- 
tion theory (hereafter SPT; d H E3, IM IH); 
a detailed description (including explicit expressions for 
P a b to third order in Pl) is presented in the Appendix. 

Comparisons with simulations (including those pre- 
sented below) have shown that the domain of applicabil- 
ity of second order (in Pj,) perturbation theory is rather 
small at z ~ 0. Furthermore, as we show below, going to 



third order is not guaranteed to improve agreement, lead- 
ing one to question the convergence properties of such a 
series expansion. If one could carry out any expansion to 
infinite order it would (trivially) give the correct answer. 
This however is usually not possible. This has led var- 
ious authors to investigate ways of summing subsets of 
the terms to arbitrary order in some expansion coefficent. 

Renormalized perturbation theory (hereafter 
RPT, see [HI, [T 3 . [2 2Tj~) is a variant of Dyson- Wyld rc- 
summation (see [1 51 ] for a discussion in the context of 
hydrodynamics) and attempts to reorganize the pertur- 
bation expansion in terms of the non-linear propagator 
and non-linear vertex to improve convergence. In partic- 
ular, if the vertex is approximated by its tree-level form 
then the power spectrum can be written as an expan- 
sion in the non-linear propagator. The resulting series is 
therefore no longer an expansion in powers of the initial 
density contrast, but rather "an expansion in orders of 
the complexity of the interaction" [23|] . 

In [Tij ] the dominant contributions to the non-linear 
propagator are identified and summed explicitly in the 
high-fc limit, giving G ab ~ e~ s k ^ 4 for large fc. Matching 
this behavior with the 1-loop propagator (valid at low k) 
gives a non-perturbative prediction for G a b- Substituting 
this propagator in the first few diagrams of the reorga- 
nized expansion then gives a non-perturbative prediction 
for the power spectrum [22j ■ We implemented the 1-loop 
and 2-loop mode-coupling contributions as described in 

The above methods work at the level of the density and 
velocity fields; an alternative approach is to use the fluid 
equations to derive equations of motion for the power 
spectrum and higher order correlators directly. Such an 
approach results in an infinite hierarchy of equations, 
which must be somehow truncated. The closure the- 
ory approach Q does so by approximating the 3-point 
correlator ($ n $t,$ r ) by its leading order expression in 
SPT. As in [lj], Gab can be computed explicitly in the 
low-fc and high-fc limits, and matched naturally in inter- 
mediate regimes. The power spectrum is then obtained 
order-by-ordcr via a Born-like series expansion. 

A variant of this approach (hereafter Time-RG the- 
ory [24|]) assumes a vanishing trispectrum to truncate 
the hierarchy. The resulting equations of motion for 
the power spectrum P a b and bispectrum B a b c can then 
be numerically integrated forward in time, starting at 
some sufficiently early redshift Zi (where P = Pl and 
B = 0). Since the time evolution is performed numer- 
ically, the method also allows the proper treatment of 
models where the linear growth factor is scale-dependent 
(e.g. models with quintessence or massive neutrinos p5j). 
This approach may be seen as a generalization of the 
renormalization group perturbation theory (here- 
after RGPT) of [2a], which is an attempt to regulate the 
relative divergence of 1-loop SPT using renormalization 
group methods. 

In [27], [H, HI] a path- integral formulation of the Vlasov 
equation is developed in terms of the distribution func- 
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tion f(x,p,t). In [3(j a similar technique is applied 
to the fluid equations (Eq. (|A40|) ). The key insight 
here is that statistical observables like the power spec- 
trum may be obtained by taking functional derivatives 
of an appropriately constructed path integral (the gen- 
erating functional). Straightforward perturbative evalu- 
ation of the generating functional reproduces the results 
of SPT, whereas applying large- N expansion techniques 
and truncating at fixed order in 1/N leads to approxima- 
tions for the power spectrum and propagator. These ap- 
proximate solutions agree with SPT up to a fixed order in 
Po, but also include non-perturbative contributions cor- 
responding to infinite partial resummations of the stan- 
dard expansion. We focus attention on the steepest- 
descent method of [3(| (hereafter Large-N), as it is con- 
siderably easier to implement than the 2PI effective ac- 
tion method. 

Lagrangian resummation theory [3l|, [32| is an ex- 
tension of the well-developed Lagrangian perturbation 
theory. Lagrangian perturbation theory (hereafter LPT; 
see [33, [3J, [351]) has received less attention recently 
than its Eulerian counterpart as a method for investi- 
gating non-linear structure growth, partly because the 
Lagrangian picture breaks down once shell-crossing oc- 
curs. However, recent work [3l( has demonstrated that 
Lagrangian perturbation theory not only reproduces the 
SPT power spectrum at the lowest non-trivial order, but 
with a slight modification also yields a non-perturbative 
prediction for the power spectrum that corresponds to 
resumming an infinite set of terms in the standard ex- 
pansion. We review LPT and the cumulant expansion in 
Appendix [B] 



III. SIMULATIONS 

In order to assess how well the perturbative expansions 
are doing, we need a reference for any given cosmology. 
We use a new set of large dynamic range N-body simu- 
lations well suited to this purpose. These computer pro- 
grams simulate the same basic physical system (a colli- 
sionless matter 'fluid' interacting only through gravity) 
that the perturbative methods attempt to describe; hence 
the results of the two methods, though arrived at very 
differently, are directly comparable. 

We have elected to investigate several different cos- 
mologies, in an attempt to better identify where and why 
various analytic techniques succeed and/or fail. For sim- 
plicity we consider only flat models in the CDM family. 
We will highlight two: the first in which a cosmological 
constant dominates the late-time evolution and which 
is close to the best-fit cosmology (ACDM: £1m = 0.25, 
fl b h 2 = 0.0224, h = 0.72, n = 0.97 and cr 8 = 0.8) and an 
extreme model (cCDM: fl M = 1, n b h 2 = 0.1, h = 0.5, 
n = 1, erg = 1) with a critical density in matter and 
a larger present-day normalization which emphasizes the 
effects of non-linearity and the erasure of baryon acoustic 
oscillations through mode coupling. 



k 




= 1) 


Ai(z = 0) 




ACDM cCDM ACDM cCDM 


0.05 


0.03 


0.03 


0.09 0.14 


0.10 


0.11 


0.09 


0.27 0.36 


0.15 


0.19 


0.22 


0.49 0.90 


0.20 


0.29 


0.27 


0.72 1.07 


0.25 


0.37 


0.38 


0.94 1.51 



TABLE I: The value of the dimensionless, linear power spec- 
trum at z = 1 and z = at several fiducial wavenumbers for 
our two example cosmologies, k is given in /iMpc -1 . 



For each cosmology the transfer function, T(k), was 
computed by evolving the coupled Boltzmann, fluid, and 
Einstein equations using the publicly available package 
CAMB (http://www.camb.info). The resulting power 
spectra were then used both as input to the perturba- 
tive methods and to generate initial conditions for the 
N-body simulations (Table H] gives the amplitude of the 
dimensionless power at some fiducial wavenumbers). 



A number of numerical issues need to be addressed in 
order to ensure that our simulations provide an adequate 
reference. Our workhorse simulations each employ 1024 3 
equal mass dark matter particles in a periodic, cubical 
box of side length 2/i _1 Gpc. By employing such large 
volumes we are highly insensitive to the periodicity of the 
box, which represents a fair sample of the Universe 0]. 
There is very little power at the fundamental mode, even 
at z = 0: A 2 (k f ,z = 0) < 10~ 4 . The lowest few modes 
obey linear growth to sub-percent accuracy and we run 
enough different realizations to ensure that the spectrum 
at the scales of interest is well determined. The large 
number of particles ensures that the spectrum is well 
converged for the fc-modes of interest, which we checked 
explicitly by comparing simulations of different box sizes. 
The simulations are evolved from Zi = 100, with the par- 
ticles perturbed from an initially uniform grid using the 
Zel'dovich approximation. The rms particle move was 
about 5% of the mean interparticle spacing. Comparison 
with second order Lagrangian perturbation theory initial 
conditions showed that this starting redshift is sufficently 
high that transients from the Zel'dovich start arc irrele- 
vant for the scales and redshifts of interest. 

Most of the evolutions were performed with a parallel 
particle-mesh code. To cross-check our results we used 
two high force resolution N-body codes: the TreePM code 
[iBl and Gadget-II [37[. These have each been tested 
against a suite of other codes @, H, Q , with very good 
agreement. We ran a subset of our simulations using all 
three codes to quantify the level of precision for the box 
size and particle loading of relevance here. With its de- 
fault time stepping, the TreePM code produces dark mat- 
ter power spectra in agreement with those from Gadget-II 
to better than 0.2% out to k ~ 1 h Mpc -1 and to O(10 -4 ) 
for k < 0.1/iMpc -1 . However these runs prove to be 
quite time consuming. If we set the time step in the 
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TrcePM code to 



2 _ 
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2 


' a ' 




0.05 


+ 


Lo.oiJ 



which evolves from 5% steps at early times to 1% steps 
as a — -> 1, we find a shortfall of power of approxi- 
mately 1% at k ~ 1/iMpc -1 but very little difference 
for fc < 0.1 /iMpc -1 . We choose the same time step- 
ping for the particle-mcsh code, which results in very 
short run times allowing an ensemble of simulations to 
be performed. With this step the particle-mesh power 
spectra show a significant deficit of power (compared to 
TreePM or Gadget-II) beyond k « 0.7/iMpc -1 but for 
k < 0.2 h Mpc -1 , the region of interest here, the agree- 
ment is better than 1%. 

To compute the power spectrum at different output 
times the particles were binned onto a regular, Cartesian 
grid using charge-in-cell assignment (38| and the result- 
ing density field was Fourier transformed. The Fourier 
modes were squared, corrected for the gridding by di- 
viding by the Fourier transform of the charge assignment 
scheme, and binned into bins equally spaced in log k. The 
average of A 2 (fc) was assigned to the average k in the bin 
and shot-noise was subtracted assuming it was Poisson. 
The binning introduces artifacts at low k, where the sam- 
pling on the grid is sparse and the dimcnsionless power 
spectrum is steep, but these are small for the scales of 
most relevance to us. Similarly there is some evidence 
that the shot-noise in simulations is not scale-invariant 
(Poisson), but the correction is negligibly small on the 
scales of interest here. 

The non-linear propagator was computed by cross- 
correlating the initial density field with the final one 
[Til ]. Similar to the power spectrum, this quantity is 
obtained by Fourier transforming both fields, multiply- 
ing their Fourier coefficients, correcting for gridding, and 
then binning the results. 

The velocity statistics are more problematic, because 
while the density and momentum fields must vanish 
where there are no tracer particles, the same is not true 
of the velocities. Thus estimates of the velocity field 
must employ a smoothing technique. Similarly the ve- 
locity field is more sensitive to finite force resolution. 
On the other hand comparison of the velocity fields with 
the density fields is less sensitive to finite volume scat- 
ter. For this reason we use a different set of simula- 
tions, with more particles (up to 3 billion) in smaller 
boxes (1.25/i -1 Gpc down to 720/i~ 1 Mpc) evolved with 
the TreePM code, for the velocity statistics. Comparison 
with different smoothing schemes, box sizes and particle 
loadings show that with these choices our results are well 
converged on the scales of interest [3!| . 



IV. COMPARISON 
A. The Power Spectrum 

We begin our analysis by comparing the predictions 
of SPT against our simulation results. Figure [T] shows 
the linear theory, 1-loop SPT, and 2-loop SPT power 
spectrum for ACDM and cCDM. While 2-loop SPT is 
a marked improvement over 1-loop SPT at z = 1, it's 
actually worse than 1-loop at z = 0. The effect is most 
apparent in cCDM, which has larger erg and This 
break-down in standard perturbation theory is not en- 
tirely surprising: since the n th order term in SPT goes 
like D 2n (z), at any given scale one expects higher order 
terms to become comparable in magnitude to lower or- 
der terms at sufficiently late times. Our results suggest 
that at BAO scales (roughly k = 0.05 - 0.25 ft, Mpc -1 ) 
the break-down occurs between z = 1 and z = 0. 

A common heuristic prescription dictates that 1-loop 
SPT can be trusted to 1% for wavenumbers satisfying 
A^(fc) < 0.4 [4l[. On the other hand a strict application 
of perturbation theory implies that 1-loop SPT can be 
trusted to 1% for wavenumbers where the 2-loop contri- 
bution is 1% of linear theory. In Figure [T] we indicate 
the predicted domain of validity of 1-loop SPT accord- 
ing to these two criteria. For comparison we also indi- 
cate where the 2-loop contribution is within 3% of linear 
theory. One sees that the agreement with simulations 
is slightly better than what our more rigorous criterion 
suggests. For instance for ACDM at z = 0, A| = 0.4 
at K w 0.12 /iMpc -1 . At this wavenumber 1-loop SPT 
overshoots the reference spectrum by about 3%, whereas 
2-loop SPT undershoots the reference spectrum by 5%. 
For cCDM at z = the situation is much worse, with 1- 
loop SPT overshooting by only 6% at fc* w 0.11 h Mpc -1 , 
but 2-loop SPT undershooting by almost 20%. 

RPT and closure theory have also been developed to 
two loops. Given the above conclusions about SPT, it is 
natural to make the same comparison between the 1-loop 
and 2-loop predictions from RPT and closure theory. In 
Figure [2] we show the matter power spectrum for these 
theories at tree, 1-loop, and 2-loop order for both ACDM 
and cCDM. For closure theory it appears that going to 
2-loop order extends the range of agreement with simu- 
lations, although the wiggles of the power spectrum are 
not matched in detail. For RPT, as with SPT, the 2-loop 
result is systematically high, whereas the 1-loop result 
performs fairly well below k ~ 0.15 /iMpc -1 . Agreement 
with simulations can be improved by changing the damp- 
ing scale in the propagator. In [22j the damping scale was 
modified by calculating S with the linear theory expres- 
sion (Eq. [7]), but using the non- linear power spectrum 
and integrating only up to k = 4fc n i. This leads to a 
~ 10% additional suppression of G(k) and hence P(k) 
on the relevant scales, brin ging the theory into better 
agreement with simulations [221 ]. At present this correc- 
tion has not been derived from first principles and we 
have not included it, but it appears that improvements 
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FIG. 1: SPT power spectrum at linear (black; dotted), 1-loop (red; solid), and 2-loop (blue; dashed) order. The squares with 
error bars show the mean and error from our N-body simulations. The four panels show ACDM (left) and cCDM (right) at 
redshifts 1 (top) and (bottom). Each curve has been divided by the no- wiggle power spectrum of [4(| to reduce the dynamic 
range. We also indicate the domain of validity of 1-loop SPT according to the heuristic prescription of [4l[ (A 2 < 0.4), and 
according to the criterion P^ < a Pl for a — 0.01, 0.03. 



in this direction could be important. 

Figure [3] shows the predicted power spectrum for the 
remainder of the theories that we consider in this work. 
With Figures [1] and [3J these figures give an overview of 
the agreement between our N-body simulations and the 
perturbation theories for ACDM and cCDM. Some of 
the trends can be seen easily in these figures, and are 
generic across cosmologies and redshifts. For instance 1- 
loop SPT, which is the same as 1-loop LPT, always over- 
predicts P(k) at high k. Lagrangian resummation theory 
on the other hand is much too strongly damped beyond 
the first wiggle. Large-./V theory more or less traces 1- 
loop SPT before turning over, while time-RG theory and 
RGPT follow the general trends of the N-body data with- 
out fitting any particular feature precisely. (Note that 
the nearly perfect agreement between RGPT and sim- 



ulations for cCDM at z = 1 is likely spurious, as this 
level of agreement is not seen for other cosmologies or at 
other redshifts.) RPT and closure give nearly identical 
tree-level predictions, and very similar 1-loop predictions 
for P(k). Closure theory appears to benefit greatly from 
going to 2-loop order, whereas for RPT even at z = 1 it 
appears that 2-loop docs worse than 1-loop. 

While we have run many realizations of each cosmol- 
ogy to reduce run-to-run variance, one sees in Figures [1] 
[5] and [3] that the N-body data are still noisy at low fc, 
which makes it difficult to make quantitative statements 
about the performance of the perturbation theories. To 
overcome this we define a 'reference spectrum' which in- 
terpolates the N-body results at high and intermediate 
k with the 1-loop SPT calculation at low fc. This elimi- 
nates the large scatter from the finite number of modes 
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FIG. 2: Comparison of the tree-level, 1-loop and 2-loop power spectrum from RPT and closure theory, for ACDM (left) and 
cCDM (right). Each curve has been divided by the no-wiggle power spectrum of [4p| ] to reduce the dynamic range. The (black) 
dotted line is linear theory, the (red) solid line is tree-level RPT, the (green) dashed line is 1-loop RPT, the (blue) long-dashed 
line is 2-loop RPT, the thick (yellow) short-long dashed line is tree-level closure, the (magenta) dot-long dashed line is 1-loop 
closure, and the (cyan) dot-dashed line is 2-loop closure. 



in the simulations and any biases from the finite bin sizes 
at low k, while still retaining the information from the 
simulations at larger k. This gives a smooth function, 
defined for all k, which can be used as a reference to 
make a quantitative comparison. Given the large num- 
ber of simulations we have run, the uncertainty in the 
N-body results is small before perturbation theory be- 
comes invalid and we can see a significant range of k for 
which theory and simulation agree well. This makes our 
final results insensitive to how the matching is done. Our 
recipe for producing a reference spectrum is to treat both 
the N-body results and 1-loop SPT as independent mea- 
surements of the true power spectrum, with errors given 



by the run-to-run variance within a wavenumber bin [54 



in the former case, and by the 2-loop SPT term in the 
latter case. Then the reference spectrum at any given k is 



defined by fitting a polynomial to all available measure- 
ments within a small wavenumber range [k — Ak, k + Afc] 
and evaluating that polynomial at k. For simplicity we 
chose to fit to a cubic with Ak = 0.01 /iMpc -1 , though 
the resulting reference spectrum is rather insensitive to 
these choices. 

All of the theories beyond linear correctly predict the 
'dip' below linear theory which can be most clearly seen 
in Figure 2] around k ~ 0.04 /iMpc -1 . This is some- 
times referred to as pre-virialization, and arises because 
the non-linear growth of the density and velocity fields is 
slower than linear on scales where the effective spectral 
index is more positive than (about) —1.5 (see jl2j for fur- 
ther discussion). In this region use of any of the methods 
provide significant improvements over linear theory 

To gain an overview of the range of validity of the 
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FIG. 3: Comparison of the power spectrum for the remaining theories. Each curve has been divided by the no-wiggle power 
spectrum of [4(j] to reduce the dynamic range. The (red) solid line is 1-loop SPT, the (magenta) dot-long dashed line is large- A 
theory, the (green) dashed line is Lagrangian resummation, the thick (yellow) short-long dashed line is time-RG theory, and 
the (cyan) dot-dashed line is RGPT. 



Methods 



Linear 

1- loop SPT 

2- loop SPT 

1- loop RPT 

2- loop RPT 

1- loop Closure 

2- loop Closure 
Time-RG 
Large-N 
Lag.Resum 



16.iZ.i8.i9, 20, 21] 
[42] 
[13, 14, 22] 

\2] 

[24] 
[30] 
[31] 



2 = 

D = 1 



0.3 
0.87 



.7 1 
72 0.63 



0.03 
0.08 
0.04 
0.10 
0.08 
0.09 
0.08 
0.04 
0.08 
0.07 



0.08 
0.10 
0.06 
0.13 0. 
0.08 0. 
0.10 0. 
0.13 
0.05 0. 
0.11 0. 
0.08 



09 0.09 
11 0.13 
09 0.23 
15 0.16 

11 0.11 
14 0.15 
17 0.27 
06 0.09 

12 0.14 
09 0.10 



various methods we list in Table [TT] the smallest value 
1-5 of k at which each method departs from our reference 
^spectrum by 1% for ACDM (a comparison with other 



H' ^sch emes defined in the literature is presented in Table 
2rjnD)- ^ s ex P ec ted, the methods perform better at smaller 
q gQScales the higher the redshift. All of the methods out- 
q ^perform linear theory, owing to the marked effects of prc- 
O.i8virialization, however none of the methods appear to be 



0.2iaccurate beyond k 
0.10 
0.17 
0.13 



0.1/iMpc" 1 at z = 0. 



B. Testing the dynamics 



TABLE II: The methods we consider in this work and the 
lowest k (in /iMpc -1 ) at which each method departs from 
our reference spectrum by 1%, as a function of redshift for 
our chosen ACDM cosmology. 



While comparison of the power spectrum is the most 
common test for perturbation theory, it is also useful 
to test if perturbation theory is correctly describing the 
underlying dynamics. To do so, we examine some of the 
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we show the cross-correlation coefficient 



0.9 



0.02 



0.04 0.06 
k (h/Mpc) 



0.08 



FIG. 4: Comparison between analytic models for P(k; z = 0) 
and the reference spectrum (Section IIV[) for model cCDM, 
focusing on large scales. Each curve has been divided by the 
no- wiggle power spectrum of to reduce the dynamic range. 
The points with error bars are the 'reference spectrum' defined 
in the text. The (black) dotted line is linear theory, the (red) 
solid line is 2-loop SPT, the (blue) long-dashed line is 2-loop 
RPT, the (green) short-dashed line is Lagrangian resumma- 
tion, the (cyan) dot-dashed line is 2-loop closure theory, the 
thick (magenta) dot-long dashed line is the large-N expan- 
sion, and the thick (yellow) short-long dashed line is time-RG 
theory. 



Ref. 


Method 


ACDM 


cCDM 






z = 


z = 1 


z = 


z = l 


[41j 


SPT 


0.12 


0.28 


0.07 


0.25 


[43] 


SPT 


0.10 


0.19 


0.08 


0.19 


® 


SPT 


0.08 


0.09 


0.08 


0.09 




Lag.R. 


0.08 


0.13 


0.07 


0.15 


® 


RPT 


0.09 


0.09 


0.09 


0.10 


[3j 


Closure 


0.09 


0.09 


0.09 


0.10 



TABLE III: The value of k, in /iMpc" 1 , to which various 
flavors of perturbation theory can be trusted according to 
various published criteria. See [|| for discussion. 



constituent pieces from the simulations, and compare to 
the perturbation theory predictions. 

Figure [5] compares the non-linear propagator 



Gi(A:) = Gii(A;) + Gi2(A;) 



(SlSI) 



(9) 



from the simulations with the predictions of analytic 
models. Only RPT and Lagrangian resummation give 
the expected behavior, G\ — > 0, for large k. 

Comparisons of perturbation theory with simulations 
typically focus on the density auto-correlation function 
or power spectrum. However perturbation theory also 
makes predictions for the (irrotational) velocity field 
which can be checked against simulations. In Figure [7] 



r(k) = Pse{k)/^Pss{k)Pee(k) 



(10) 



for several theories, compared with the same quantity 
measured from simulations. In linear theory r{k) = 1 
identically. On physical grounds one expects to see a de- 
coherence of density and velocity fields on small scales, 
and indeed the simulations show r(k) — > for large k. 
None of the analytic theories correctly reproduce this be- 
havior. SPT and Time-RG theory follow the downward 
turn of the simulation data initially, but then predict an 
unphysical r(k) > 1 very soon after non-linear correc- 
tions become important. RPT and closure theory per- 
form somewhat better, in that r(k) never exceeds unity, 
but the level of agreement with simulations is still not 
good above k ~ 0.1 h Mpc -1 . (Note that we have dis- 
played here only the 1-loop predictions from these theo- 
ries.) The deviation in r(k) seems to be driven mostly by 
the densities, with perturbation theory performing bet- 
ter at the same scale for the velocities than the densities 
(see Figure E]). 



V. DISCUSSION 

Standard perturbation theory has a simple and direct 
theoretical motivation, and results in explicit integral ex- 
pressions at any order. If taken to infinite order, it pro- 
vides an exact solution (though to an idealized problem). 
While standard perturbation theory works well at high 
redshift and large scales, our results indicate that the 
standard expansion is badly behaved at the redshifts and 
scales most accessible to observation, in that higher or- 
der terms are comparable in magnitude to lower order 
terms. Although one expects the expansion to converge 
if taken to sufficiently high order, this comes at a great 
computational cost. With advances in raw computing 
power it may one day become possible to perform the 
calculation to the requisite order, but in the near future 
this approach seems impracticable. 

On the other hand, it should be emphasized that SPT 
performs rather well at high redshifts, z > 1. Figure 
[1] shows that 2-loop SPT at z = 1 agrees with simula- 
tions to 1% out to k = 0.2/iMpc -1 or beyond (where 
the simulations themselves become unreliable). At these 
redshifts SPT not only provides a reasonable theoretical 
prediction for the matter power spectrum on observation- 
ally relevant scales, but also an estimate of the theoretical 
uncertainty on this prediction. 

RPT is essentially a rearrangement of the standard 
expansion, so like SPT it is an exact solution if car- 
ried out to all orders. While this rearrangement appears 
to improve the convergence properties of the perturba- 
tion series, it makes it unclear what small quantity (if 
any) we are actually expanding in. Furthermore, RPT 
does not actually provide closed-form expressions for the 
power spectrum, but rather integral relations where P a b 
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0.15 



FIG. 5: The fractional deviation of each method from the reference spectrum, for ACDM at z — (left) and z = 1 (right). 
This figure focuses on the region k < 0.15 /iMpc -1 where linear theory is inadequate but higher order methods are still viable. 
As in Figure |4] the (black) dotted line is linear theory, the (red) solid line is 2-loop SPT, the (blue) long-dashed line is 2- loop 
RPT, the (green) short-dashed line is Lagrangian resummation, the thick (cyan) dot-short dashed line is 2-loop closure theory 
the thick (magenta) dot-long dashed line is the large-A expansion, and the thick (yellow) short-long dashed line is time-RG 
theory. 




0.1 0.2 0.3 0.4 



k (h/Mpc) 

FIG. 6: The non-linear propagator (normalized to 1 at k = 0) 
for ACDM at z = 0. The (red) solid line is SPT, the (green) 
short-dashed line is Lagrangian resummation, the (blue) long- 
dashed line is RPT, the thick (cyan) dot-short dashed line is 
closure theory and the thick (magenta) dot-long dashed line 
is the large- A expansion. 



is expressed in terms of mode-coupling integrals of it- 
self. Thus in addition to truncating the loop expansion 
at finite order, a fully consistent implementation of RPT 
requires solving for P a b according to an iterative scheme, 
of which the explicit expressions presented in [2^ | repre- 
sent only the first step. The error associated with this 
approximation has (to our knowledge) yet to be quanti- 



fied. 

Closure theory derives from a very different pcrturba- 
tive scheme than RPT, yet the results obtained are super- 
ficially quite similar. There is no obvious way to provide 
error estimates on the results of closure theory, however, 
as the closure equations are obtained from heuristic ap- 
proximations rather than a systematic expansion. Fur- 
thermore the propagator in closure theory shows unre- 
alistic oscillations for large k. As mentioned previously, 
the closure equations are solved approximately in 0] by 
means of a Born-like expansion. Recently [44j an attempt 
has been made to solve the closure equations numerically 
without resort to such a Born-like expansion. The result- 
ing predictions for the power spectrum appear to agree 
better with simulations than the results presented here, 
although it is difficult to draw any firm conclusions from 
the information provided. 

Time-RG theory is based on a single well-defined ap- 
proximation: the vanishing of the trispectrum. The va- 
lidity of this approximation can easily be checked in sim- 
ulations, and in principle this could allow one to quantify 
the theoretical uncertainty in the method. As most eas- 
ily seen in Figure although time-RG theory follows 
the general trends of our reference spectrum over a wider 
range than other methods, it comes up short by 1-2% over 
the entire quasi-linear regime. It also gives an unphysical 
prediction for the density-velocity cross-correlation. 

The large- N expansion utilizes more sophisticated the- 
oretical machinery than other resummation techniques. 
While the path-integral formalism for computing cluster- 
ing statistics is exact, the errors introduced by the largc- 
N expansion are difficult to quantify, as W is a fictitious 
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FIG. 7: A comparison of the density-velocity cross-correlation 
predicted analytically with that measured in simulations, for 
ACDM at redshift z = (top) and z « 1 (bottom). As 
in Figure [4] the solid (red) line is SPT, the dashed (blue) 
line is RPT, the dot-dashed (cyan) line is closure theory, and 
the short-long-dashed (yellow) line is Time-RG theory. For 
simplicity we show only the 1-loop predictions for SPT, RPT, 
and closure theory. 



parameter. Although the large- TV expansion corresponds 
to an infinite partial resummation of the standard per- 
turbativc expansion, from our results it seems that this 
resummation offers little improvement over 1-loop SPT 
in the quasi-linear regime. The grossly unphysical behav- 
ior of the propagator in this theory is likely responsible 
for this effect. As mentioned previously, we have focused 
attention on the steepest-descent method rather than the 
2PI effective action method of [3(J. The latter method 
produces a more reasonable propagator, and likely results 
in a better prediction for the power spectrum, although 
at an increased computational cost. 

Like SPT, the Lagrangian resummation prescription 
of [3l| also results in easy to compute, explicit integral 
expressions. These are well behaved at large k, allowing 
e.g. £(r) to be computed, and there are natural extensions 
to redshift space and to halo bias [32j . For the real-space 
mass power spectrum considered here it offers a marginal 
improvement over 1-loop SPT for fcS < 1/2, although 
the damping prefactor strongly overcompcnsates as one 
moves further into the quasi-linear regime. 

Our results have interesting implications for generating 
a suite of simulations aimed at constraining the matter 
power spectrum. If we can trust perturbative methods 
for fc£ < x c , then we can focus the computational re- 
sources on higher k. Assuming Gaussian fluctuations, 
obtaining 1% accuracy in a bin (fc; Afc) requires 2 x 10 4 



k 


A? in (fc) 


A r 2 cf (fc) 


G(k) 


0.02 


0.012 


0.012 


0.996 


0.04 


0.053 


0.053 


0.980 


0.06 


0.130 


0.129 


0.950 


0.08 


0.210 


0.210 


0.914 


0.10 


0.274 


0.285 


0.859 


0.12 


0.398 


0.410 


0.804 


0.14 


0.466 


0.507 


0.737 


0.16 


0.533 


0.617 


0.664 


0.18 


0.662 


0.764 


0.592 


0.20 


0.720 


0.894 


0.518 



TABLE IV: Our input linear theory spectrum, at z = 0, for 
the ACDM model as a function of wavenumber (in ftMpc -1 ) 
and the reference spectrum and propagator [G(fc)] from our N- 
body simulations. Our convergence tests indicate the spectra 
should be accurate to < 1% over the range of scales shown. 

modes. There are (fcLbox) 3 (Afc / fc) / (2tt 2 ) modes from a 
periodic box of side length Lbox, so our 1% constraint at 
fc£ ~ x c translates into 

S f 2n 2 N \ 1/3 
ibox ~ x c \ Ak/k ) 

* 3Gpc © (lost) (^y ( Ai7fc) n) 

or an equivalent volume of smaller simulations. This con- 
straint is most difficult to meet at z = 0, since £ is larger 
and the simulations must be evolved for longer. As an ex- 
ample with the default parameters listed above we would 
require 27 simulations, each 1 h~ 1 Gpc on a side, to ob- 
tain percent level constraints on the power spectrum of 
ACDM in a 10% band near k ~ 0.1 h Mpc" 1 at z = but 
at z = 1 we could trust perturbation theory at this scale 
and focus the simulations on k ~ 0.15 h Mpc -1 where 
three times fewer simulations of the same size are needed. 



VI. CONCLUSIONS 

Perturbative methods have a long history in cosmol- 
ogy, and are widely used in many fields of physics. Many 
of the techniques reviewed herein were first developed in 
other fields and applied to other problems, with vary- 
ing levels of success, before being pressed into service for 
modeling cosmological perturbations. In this paper we 
have studied a variety of these methods as applied to 
predicting the large-scale clustering of cold, collisionless 
matter in an expanding Universe. Our results indicate 
that the analytic theories correctly model the approach 
to non-linearity and work well when the perturbations 
are small, but none of the available theories are up to 
the challenge of fully describing the behaviour of matter 
on quasi-linear scales at late times. We have emphasized 
the need to study a range of different cosmologies and 
to look at a variety of different statistical observables, as 
accidental agreement between theory and simulations is 
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FIG. 8: The density-velocity cross spectrum (left) and the velocity-velocity auto-spectrum (right) for the ACDM model at 
z = 0. As in Figure [7] the (black) dotted line is linear theory, the (red) solid line is 1-loop SPT, the (blue) long-dashed line 
is 1-loop RPT, the (cyan) dot-short dashed line is 1-loop closure theory, the thick (magenta) dot-long dashed line is large- N 
theory, and the thick (yellow) short-long dashed line is time-RG theory. 



possible if one only considers the power spectrum. We 
have computed the 2-loop contribution to SPT and found 
that the standard pcrturbativc expansion is badly be- 
haved at low redshifts, even on scales where 1-loop SPT 
was previously believed to be valid. This provides further 
motivation for studying alternative analytic approaches 
based on non-perturbative methods, though at the same 
time it emphasizes the need for error control in analytic 
methods. 

This work has made use of a large number of high dy- 
namic range N-body simulations, against which we can 
compare the analytic models. We make these data public 
in Table HVl to aid future work in the field. In addition 
a flexible software package that implements the pertur- 
bation schemes described in this paper is available from 
the authors. 
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APPENDIX A: EULERIAN PERTURBATION THEORY 



Here we briefly recap the derivation of the fluid equations in the Eulerian picture, and the assumptions that are 
made in perturbative treatments ( [H, 03, 0, EH HO, Hl[ ; see [l2| for a review). The matter content of the Universe is 
modeled as a large collection of identical particles of mass m, interacting only through mutual gravitational attraction. 
For low densities and sub-horizon scales, such forces are adequately described by Newtonian gravity in a uniformly 
expanding background, with the Newtonian potential sourced by inhomogeneitics in the density field. The distribution 
function for such a set of particles obeys the Vlasov equation. The N-body methods are essentially a Monte-Carlo 
evolution of the Vlasov equation where the Monte-Carlo tracer super-particles move along characteristics. 

Analytically one typically invokes the single-stream approximation, which assumes that all particles at a given point 
x move together with the same velocity v{x). This amounts to demanding that f(x,p) oc Sr>\p — mav(x)], where / is 
the distribution function and 5d is the Dirac delta function. This assumption is explicitly violated once shell crossing 
occurs in gravitational collapse, but is thought to be a reasonable approximation for small density constrasts. The 
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velocity moments of the Vlasov equation then give the familiar fluid equations (e.g. 

^ + V-[(l + 6)v]=0, (Al) 
-^-+Hv + iv- V)v + V$ = 0. (A2) 

OT 

where ri = dhia/dr = aH is the conformal Hubble parameter. 

It is conventional to further assume that the vorticity w = V x v of the velocity field vanishes, i.e. that the fluid 
is irrotational. This assumption is motivated by noting that w oc a" 1 at linear order, and is well supported by 
simulations |45l |46|. Under this approximation the velocity field is completely specified by its divergence 9 = V ■ v, 
and the fluid equations reduce to 

!*+0 = -V-(<fo), (A3) 

OT 

r)f) 

— + HO + 4nGa 2 p5 = - V • [(« • V)v]. (A4) 
In Fourier space v(k) = —ikd(k)/k 2 , giving 

« + neik) + ^ m H 2 <5(fc) = -/ ^^fofe + 92 - (A6) 

As long as (5 and are small, the right-hand sides of the fluid equations are small and can be dropped; this 
approximation defines linear theory. The solution to the resulting linearized fluid equations may be written as 

S L (k;z) = ^\s i (k) , e L (k;z) = -H(z)f(z)^\s i (k), (A7) 
D{zi) D{zi) 

where Si is the density contrast at some early time Zj when linear theory is certainly valid, D is the linear growth 
function (normalized to 1 today), and f = dhiD/dliia. At early times Q m rj 1 and D oc a. Note that a possible 
decaying mode contribution, proportional to a" 3 / 2 at early times, is forced to zero in linear theory by the condition 
that 5 be well-behaved as a — > 0. Note also that the mode-coupling integrals vanish for k = 0, so linear theory is 
always valid in some neighborhood of k = 0, even at late times. For convenience we define to be the linear density 
contrast today, i.e. So(k) = Si(k;z = 0). 

It often proves convenient to use rj — lnl? as a time variable, and to combine <5 and into a two-component field 

*-(*)=(_ JSLt)- (A8) 



-0(k)/Hf 



If we introduce 



kqi a( . k 2 (q 1 -q 2 ) 



a(9i,92) = — — > P(9i>92) = — ^xi — ( A9 ) 



the fluid equations may be recast as 

x 9 o 
a?7 



$ 6 (fe; 7?) = y ^3 7 ahc (g, fc - g)* 6 (g; ?y)$ c (fc - 9; »?), (A10) 



where 



Haifa) = ( _3 n ra SlU 1 .! ) (All) 



2F 2p 

and the vertex 7 & c (gi,92) only has nonzero entries 7121(91,92) = 7112(92,91) = a(qi,q 2 )/2 and 7222(91,92) = 
/3(<7i, <72). The initial fields at time 77, are denoted 

0„(fc) = * o (fc;»fc) = *i(*o(j), (A12) 



and the linear theory solution is simply $ a ; (fe; 77) = e n 71i (j) a {k) 
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FIG. 9: Diagrammatic representation of the n th order contribution to S(k). 



1. Beyond linear order 

Standard perturbation theory (hereafter SPT; [l6|, [13, EH EM H3, HH) defines a systematic series solution to the 
fluid equations ([T]) in powers of the initial density contrast Si (or equivalently in powers of the current linearly evolved 
density contrast 5q). In an Einstein-de Sitter universe, where TL oc a -1 / 2 and Sl m TL 2 oc a" 1 , the expansion may be 
written as 

oo oo 

5(fe;r) = ^a™(r)5„(fe), fl(fc;r) = -W(r) £o n (r)fl„(fc), (A13) 

n— 1 n— 1 



where S n (k) and n (k) are time-indepedent mode-coupling integrals over n powers of the initial density field: 

■ / %i^< 2 *>^ d> - *) (5!!) *(«) ■ ■•*-(*)■ ( a ») 



The kernels -F n and G ra satisfy recurrence relations that follow straightforwardly from the equations of motion [171 . 
EE HI: 



j?, _ ^4 G m {q\, . . . , q m ) r fe ■ fci fc 2 (fci-fc 2 ) . 

^nl/Zl) ■ • ■ , 9nJ — 7. TTJ TV U + Zn ) U 2 ^n-m\,Qm+l, ■ ■ ■ , Qn) H TJT2 ■ ■ • , <ZnJ 

(A15) 

^ , v G m (qi, . . . ,q m ) r fe ■ fei fc 2 (fei-fe 2 )„ , >] , 

<~*n\qi, ■ ■ ■ ,q n ) — 7^ — r^jT? tt j_ FF J '»-" i W m + 1 '-"'W + T1 — rrrs — ^n-mi?m+i,---,9nj , l Ai oj 

m=1 (zn + 6)(n - L) i k 1 K 1 K 2 J 

where k± = qi + h q m , k 2 = q m +i H h q„, fc = fei + fc 2 and F\ = G\ = 1. 

While the Einstein-de Sitter approximation is convenient, it is not necessary |47j |. However we have confirmed that 
an accurate approximation is to substitute the growth factor D(z) for a, 

oo oo 

6(k;z) = 2Z Dn ( z ) S n(k), 9(k;z) = -H.f2Z Dn ( z ) d n(k), (A17) 

n— 1 n— 1 

with the same mode-coupling integrals as above for S n and 6 n . The validity of this approximation is ultimately traced 
to the fact that the ratio fl m / f 2 is very nearly unity over the entire lifetime of the universe for ACDM cosmologies, 
since f^Q° m 6 M- 

To compute statistical observables it is convenient to introduce diagrammatic rules for keeping track of the various 
terms in the perturbation series [F7| . The function 5 n (k) (or 9 n (k)) may be represented as in Figure^ where the open 
circles denote factors of 5q, and the vertex denotes a momentum-conserving integral of F n (or G n ) over intermediate 
wavevectors q^. Algebraically the n th order contribution PW is obtained by isolating all terms of order (Jo) 2 " from 
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FIG. 10: Diagrammatic prescription for computing P*- 2 ' 2 ' (£:). The overall factor of 2 comes from the two equivalent ways 
of pairing the open circles. Only the single wavevector q must be integrated over, the rest being determined by momentum 
conservation at vertices and translational invariance of the 2-point function. 



the ensemble average (8(k)8(k')), i.e. 

2n-l 

(2n) 3 5 D (k + k')P^(k; z) = D 2n (z) (S m (k)d2n- m (k')). 



(A18) 



The quantity {S m (k)S2n-m(k')) may be represented diagrammatically by "multiplying" the diagrams for S m (k) and 
S2n-m(k'). Since the initial field Si (and hence do) is Gaussian, ensemble averages of powers of So may be expanded in 
terms of the 2-point function Pq according to Wick's theorem. Then the product of the diagrams S m (k) and S2 n -m(k') 
is given by summing over all possible pairings of their open circles, where open circles are paired according to the rule 



-O x o- 



(2w) 3 6 D (q + q'Wq), 



(A19) 



with the additional understanding that any diagram containing a tadpole (a fragment connected to the rest of the 
diagram by a single edge) vanishes identically. 

As an example we show in Figure [TO] how to obtain the 2nd order contribution p( 2 ' 2 )(fc). Notice that after invoking 
momentum conservation at vertices and translational invariance of the 2-point function, only a single wavevector 
remains to be integrated. In general all diagrams contributing to p(") contain n — 1 loops, requiring integration over 
n — 1 independent wavevectors. For this reason we often classify power spectrum terms by their number of loops 
rather than their "order," which is a potentially ambiguous concept. 

With this expansion, statistical observables may be computed straightforwardly in SPT to any fixed order. For 
example, the first correction to the matter power spectrum (second order in the initial power spectrum, fourth order 
in the initial density contrast, or 1-loop in the diagrammatic idiom) is given by 



P(k) = P L (k) + P< 2 ' 3 )(fc) + P^(k) 



where Pt(fc; z) = D 2 {z)Po(k) is the linear power spectrum and pi 



(kr) 
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2rxy 



(A20) 

(A21a) 
(A21b) 



At low k, p( 2 ' 2 ) is positive while P^ 1 ' 3 ) is negative, and there is a large degree of cancellation between them. For large 
k 



pW(k)~^Y?k 2 P L {k) and P^ 3 \k) ~ -i£ 2 fc 2 P L (fc) 

where S is defined by Eq. ([7]), so for sufficiently large k the total second order contribution is negative. 
It is also straightforward to derive expressions for the velocity power spectrum [l8j 
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(A22) 

(A23) 
(A24) 
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the density-velocity cross-spectrum (e.g. [4£ 



P^ 3 \k) = ~P L (k) / dr P L 

59 V ' 252 4tt 2 W io 



24 3 
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98 4^ 
and the propagator 



dx P L (k\/l + r 2 - 2rx 
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(1 



2rxf 



G{k) ~ 1 + 



p(l,3) + p(l,5) 

2Pl 



(A25) 
(A26) 

(A27) 



Though it is not usually considered, there is no real obstacle in going to the next order in the systematic perturbative 
expansion described above. For the third order (2-loop) contribution one finds P^ 3 '(k) = P^' 5 )(fc) + P^ 2 ^(k) + 
p( 3 ' 3 )(fc) with [H 



pW(k)=30P L (k) J F^(k,q,-q,p,-p)P L (q)P L 



(P) 



V ' J (2tt) 3 (2tt)3 



(2tt) 3 (2tt) 3 

F ( 2 s \ qi k- q )Fl s) {-q,q- k,p, -p)P L (q)P L (p)P L (\k - q\) 
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(A28) 
(A29) 

(A30) 



and with given by Eq. (|A15|) symmetrized over its n arguments qi, . . . ,q n . Using rotational symmetry to 
eliminate one azimuthal integration, the resulting expressions require 5-dimensional mode-coupling integrals which 
are best performed using Monte Carlo methods. 

Expressions for higher order contributions are not difficult to derive, but the computational costs of evaluating 
them quickly spiral out of control. In general the £-looj> contribution requires mode-coupling integrals of dimension 
3£ (3£ — 1 after rotational symmetry), making 1-loop simple, 2-loop possible, and higher orders impracticable. 



APPENDIX B: LAGRANGIAN PERTURBATION THEORY 



The Lagrangian description of structure formation [49l. l50l. l5l| relates the current, or Eulerian, position of a mass 
element, x, to its initial, or Lagrangian, position, q, through a displacement vector field: x = q + ^f(q). (Note that 
q is used as a position vector in the Lagrangian picture, whereas the same symbol is used as a wavevector in the 
Eulerian picture.) The displacements can be related to overdensities by [52| 



5{x) = J d 3 q5 D {x-q-V)-l , 5(k) = J d 3 q e- lkq (e~ ik *W - 1 
The displacements evolve according to 



d 2 * 

~dt 2 ~ 



d^f 

2H— = -V a 
dt 



*(«?)] 



(Bl) 



(B2) 



where here and only here <p is the gravitational potential. Analogousto Eulerian perturbation theory, standard LPT 
expands the displacement in powers of the linear density field with [53| 



n 

i=l 



d 3 h 



(2k) 3 5 D \J2 k i- k ) i(n) (fei, • • • ,kn, k)6o(ki) ■ ■ ■ 5 (k n ) 



(B3) 



and the lS n > have closed form expressions in terms of dot products of wave vectors which can be generated by 
recurrence relations. Expanding the exponential in Eq. (|F31[) we obtain a perturbative series for the overdensity, 
6 = <S (1) + <5 (2) H where, e.g., 



^ 2 )(fc) 



d 3 k t d 3 k 2 
(2tt) 3 



S D (ki +k 2 - fc)5 (fci)5o(fc 2 ) \k • L (2) (fci,fc 2 , k) + k ■ ■ L^{k 2 )] (B4) 
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is second order in the linear density field So- 

A similar expansion can be performed for the power spectrum, which from Eq. (jBip can be written 



P(k) 



d 3 q e~ ik q 



« 



— ifc-A* 



(B5) 



where A* = VP (q) — *(0) and we have used translational invariance. 

Alternatively [3l| suggested using the cumulant expansion theorem for the exponential in Eq. (|B5|) and using the 
binomial theorem to expand the term (k ■ A*&) N . One obtains two types of terms: those depending on \& at the same 
point and those depending on \& at two different points. Owing to statistical homogeneity the first type of term is 
independent of position and can be factored out of the integral leaving [3l[ 



P(k) = exp 



exp 



^ Nl 

JV=2 



- 1 



(B6) 



where fc,-. 



k- B [N) 



(r) is shorthand for the second type of term. 



In a traditional pcrturbative calculation one would expand this expression to a fixed order in \&; this approach 
indeed reproduces the SPT result to 2nd order. However one might expect that the position-independent cumulant 
factors are more important on large scales than the position-dependent ones, suggesting that these factors should be 
left unexpanded in the exponential. Using well-known previous results from LPT the first corrections to the power 
spectrum are then [3l[ 



P(k) = e -( fcS ) 2 /2 rp L (fc) + p(2.2)(jfe) + P(b3)( fc ) 

where S is given by Eq. 0, p( 2 < 2 >(fc) is as in SPT [Eq. (|A21b|) ] and 



{kr) 



if + 10 + 100r 2 - 42r 4 + -^-(r 2 - l) 3 (7r 2 + 2) In 



(B7) 



(B8) 



Notice that this differs from the SPT result [Eq. (|A21a[) ] only by the replacement —158 — > 10 in the brackets. If the 
exponential prefactor is expanded to first order in Pj, the SPT result is recovered exactly [3l[ . Also note that the first 
term k ' 2 PT,(k) is identical to the tree-level RPT result in the large-fc limit. 
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